# Read in Data

 all_mfri = list.files('./Wildfire_MFRI/',pattern = '.tif',full.names = T)
 for(file in all_mfri){
   object_name = file_path_sans_ext(basename(file))
   assign(object_name, raster(file))
 }

CC4a_reg = read_sf('./Boundries/CC4a_RegionsSub.shp')
CC4a_reg  = st_transform(CC4a_reg, "+proj=aea +lat_1=34 +lat_2=40.5 +lat_0=0 +lon_0=-120 +x_0=0 +y_0=-4000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
# dissolve multipart Sierra Mountain feature
CC4a_reg = ms_dissolve(CC4a_reg,field = 'Region')

# Get stats for 01-25 & 26-50 MFRIs cap at 500 years

 all_2001_2025 = stack(all_mfri[grepl('2001_2025',all_mfri)])
 all_2026_2050 = stack(all_mfri[grepl('2026_2050',all_mfri)])
 
 # years for maximum plotted MFRI
 capat=500

 # calculate stats for stacks and cap MFRI 
 # writes out files with 3 different postfix _stat for statisic (mean etc), _df for dataframe for ggplot, _capped for raster capted at value capat
 summary_functions = c('min','max','mean' )
 for(summary in summary_functions){
   for(stac in c('all_2001_2025','all_2026_2050')){
      assign(paste(summary,stac,sep='_'),do.call(summary,list(x=get(stac),na.rm=T)))
      capper = get(paste(summary,stac,sep='_'))
      # write out df for ggplot
      capper.df = data.frame(rasterToPoints(capper))
      names(capper.df)  =c("lon", "lat","MFRI")
      assign(paste(summary,stac,'df',sep='_') ,capper.df)
      # cap at catat yrs
      capper[capper>capat]=capat 
      assign(paste(summary,stac,'capped',sep='_'),capper)
   }}
 
  # write out df for ggplot
  MFRI_76_00.df = data.frame(rasterToPoints(MFRI_76_00))
  names(MFRI_76_00.df)  =c("lon", "lat","MFRI")
all_2001_2025[all_2001_2025>500]=500
all_2026_2050[all_2026_2050>500]=500
plot(all_2001_2025)

plot(all_2026_2050)

 mean_chg_76_25 = MFRI_76_00.df
 mean_chg_76_25$MFRI =  mean_all_2001_2025_df$MFRI-mean_chg_76_25$MFRI
 mean_chg_76_25$MFRI[ mean_chg_76_25$MFRI >350]=350
 rng= range(mean_chg_76_25$MFRI)
 
 ggplot()+geom_raster(data=mean_chg_76_25,aes(x=lon,y=lat,fill=MFRI))+ 
   scale_fill_gradientn(colours= c("#cc0000", "#cc0000"  , 'grey', "#339933","#339933" ), #colors in 
      limits=c(-350, 350))+ geom_sf(data=CC4a_reg,colour = "grey30", fill = NA,size=.75) + #same limits for plots
   ggtitle('Change in MFRIs 2000 - 2025 mean model run')+coord_sf()

for(region_aoi in CC4a_reg$Region){  
    aoi = st_bbox(CC4a_reg[CC4a_reg$Region==region_aoi,])
    aplot = ggplot()+  geom_raster(data=mean_chg_76_25,aes(x=lon,y=lat,fill=MFRI))+ 
     scale_fill_gradientn(colours= c("#cc0000", "#cc0000"  , 'grey', "#339933","#339933" ),  
        limits=c(-350, 350))+ geom_sf(data=CC4a_reg,colour = "grey30", fill = NA,size=.75) +  
     ggtitle(paste(region_aoi,'\nChange in MFRIs 2000 - 2025 \nMean model run'))+
      coord_sf( xlim=c(aoi[1],aoi[3]),ylim=c(aoi[2],aoi[4]) )
    plot(aplot)
    }

 mean_chg_76_50 = MFRI_76_00.df
 mean_chg_76_50$MFRI =  mean_all_2026_2050_df$MFRI-mean_chg_76_50$MFRI
 mean_chg_76_50$MFRI[ mean_chg_76_50$MFRI >350]=350
 rng= range(mean_chg_76_50$MFRI)
 
 ggplot()+geom_raster(data=mean_chg_76_50,aes(x=lon,y=lat,fill=MFRI))+ 
   scale_fill_gradientn(colours= c("#cc0000", "#cc0000"  , 'grey', "#339933","#339933" ),   
                 limits=c(-350, 350))+ geom_sf(data=CC4a_reg,colour = "grey30", fill = NA,size=.75)+  
    ggtitle('Change in MFRIs 2000 - 2050 mean model run')+coord_sf()

for(region_aoi in CC4a_reg$Region){  
    aoi = st_bbox(CC4a_reg[CC4a_reg$Region==region_aoi,])
    aplot = ggplot()+  geom_raster(data=mean_chg_76_50,aes(x=lon,y=lat,fill=MFRI))+ 
     scale_fill_gradientn(colours= c("#cc0000", "#cc0000"  , 'grey', "#339933","#339933" ),  
        limits=c(-350, 350))+ geom_sf(data=CC4a_reg,colour = "grey30", fill = NA,size=.75) +  
     ggtitle(paste(region_aoi,'\nChange in MFRIs 2000 - 2050 \nMean model run'))+
      coord_sf( xlim=c(aoi[1],aoi[3]),ylim=c(aoi[2],aoi[4]) )
    plot(aplot)
    }

 the_fun = median  # function used to summarize raster values by polygons
 region_code = data.frame(ID=seq(1,9),as(CC4a_reg,'Spatial')@data$Region)
 full_mean_stack = stack(MFRI_76_00,mean_all_2001_2025,mean_all_2026_2050)
 names(full_mean_stack) =c('2000','2025','2050')
 extract_full_mean_df = extract(full_mean_stack, as(CC4a_reg,'Spatial'), fun=the_fun, na.rm=T, df=T)
 extract_full_mean_df = left_join(region_code,extract_full_mean_df,by='ID') %>% select(-ID)%>%melt()
 extract_full_mean_df$Year = as.numeric(substr(as.character(extract_full_mean_df$variable),2,5))
 names(extract_full_mean_df)=c('Region','variable','value','Year')

 the_fun = min  # function used to summarize raster values by polygons
 full_min_stack = stack(MFRI_76_00,min_all_2001_2025,min_all_2026_2050)
 names(full_min_stack) =c('2000','2025','2050')
 extract_full_min_df = extract(full_min_stack, as(CC4a_reg,'Spatial'), fun=the_fun, na.rm=T, df=T)
 extract_full_min_df = left_join(region_code,extract_full_min_df,by='ID') %>% select(-ID)%>%melt()
 extract_full_min_df$Year = as.numeric(substr(as.character(extract_full_min_df$variable),2,5))
 names(extract_full_min_df)=c('Region','variable','value','Year')
  ggplot()+geom_smooth(data=subset(extract_full_min_df, Region != 'Inland South' ),aes(x=Year,y=value,color=Region,group=Region),size=1.25,alpha=.6)  +ggtitle('Minimum observed MFRI of Min(all_models) by region - \n omitting Inland South')

  ggplot()+geom_smooth(data=subset(extract_full_mean_df, Region != 'Inland South' ),aes(x=Year,y=value,color=Region,group=Region),size=1.25,alpha=.6)  +ggtitle('Median observed MFRI of Mean(all_models) by region - \n omitting Inland South')